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ABSTRACT 

The steady-state solution of the system of equations consisting of the full Navier-Stokes equa- 
tions and two turbulence equations has been obtained using a multigrid strategy on unstruc- 
tured meshes. The flow equations and turbulence equations are solved in a loosely coupled 
manner. The flow equations are advanced in time using a multi-stage Runge-Kutta time step- 
ping scheme with a stability bound local time-step, while the turbulence equations are 
advanced in a point-implicit scheme with a time-step which guarantees stability and positivity. 
Low Reynolds number modifications to the original two-equation model are incorporated in a 
manner which results in well behaved equations for arbitrarily small wall distances. A variety 
of aerodynamic flows are solved for, initializing all quantities with uniform freestream values. 
Rapid and uniform convergence rates for the flow and turbulence equations are observed. 
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1. INTRODUCTION 

The use of unstructured meshes has become more widespread in recent years due to the 
ease with which complex geometries can be handled and the possibility of enhancing the solu- 
tion accuracy and efficiency through adaptive meshing techniques. To date, most of the 
successes of unstructured mesh techniques have been in computing inviscid flows in two and 
three dimensions over arbitrary geometries. However, more recently, solutions of the Navier- 
Stokes equations on unstructured meshes have been reported [1,2, 3, 4, 5]. The main obstacles to 
efficiently computing high-Reynolds-number flows on unstructured meshes are due to the 
required grid stretching and the turbulence model. For high-Reynolds-number flows over 
streamlined bodies, viscous effects are confined to thin boundary-layer and wake regions, 
which can only be resolved efficiently using high aspect ratio elements. One approach [3,5] is 
to fit a thin local mesh of structured high aspect ratio quadrilaterals in the viscous regions, and 
fill the remainder of the domain with an unstructured mesh. The other approach consists of 
filling the entire domain with an unstructured mesh which contains highly stretched triangular 
elements in the viscous regions [4], In this work, the latter approach has been pursued, in the 
interest of developing a more general method capable of dealing with a wider variety of flows, 
such as flows with confluent boundary layers, or mixing wakes, and also to enable the 
straight-forward implementation of adaptive meshing techniques throughout all regions of the 
flow-field. The numerical scheme must therefore be formulated such that the accuracy and 
convergence are not seriously affected by the presence of highly stretched triangular elements. 

The most commonly employed turbulence models for compressible flow calculations are 
of the algebraic mixing-length type [6]. These models have been shown to produce good 
results for attached turbulent boundary layers and mildly separated flows using structured 
meshes, and have also been implemented for non-trivial geometries on unstructured meshes [7]. 
Although such models can be made inexpensive and computationally robust even in the context 
of unstructured meshes, they lack the generality required for dealing with completely arbitrary 
geometries, and their ability in predicting flows with multiple confluent shear layers and large 
amounts of separation is at best limited. Two equation models, on the other hand, offer the 
possibility of dealing with the more complicated flows which are often associated with the 
complex geometries for which unstructured meshes are so well suited. In principle, the imple- 
mentation of such models on unstructured meshes can be accomplished in a straight-forward 
fashion, simply by discretizing and integrating the turbulence equations in a manner analogous 
to that employed for the mean flow equations. However, field-equation turbulence models have 
often proved to be extremely difficult to integrate to steady-state, exhibiting stiff or unstable 
numerical behavior in regions very close to the wall, as well as in the far-field. The use of 
multigrid to solve the turbulence equations has recently been reported by several authors [8,9], 
using a Ni-type scheme on structured meshes. In this work, a multigrid strategy which has 
previously been developed for the Euler and Navier-Stokes equations on unstructured meshes 
[4,10] is extended to solve for the two turbulence equations as well. 

2. GOVERNING EQUATIONS 

The governing equations are obtained by Favre averaging the Navier-Stokes equations, 
and modeling the Reynolds stress and heat flux terms by the Boussinesq assumption. In con- 
servative form, these equations are written as 
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where w is the solution vector and f c and g c are the cartesian components of the convective 
fluxes 
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In the above equations, p represents the fluid density, u and v the x and y components of fluid 
velocity, E the total energy, and p is the pressure which can be calculated from the equation of 
state of a perfect gas 


P = (Y-l)p 

The viscous fluxes /„ and g v are given by 


E - 


(i^+v 2 ) 


(3) 



'o 


'o 



8v = 




°W 


K< 


. u °>* +v °yy . 


(4) 


where a represents the stress tensor, and q the heat flux vector, which are given by 


= 2(p + \i,)u x - j (p + p,)(u,+v,) - -|p* 
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P represents the molecular viscosity, and p, denotes the turbulent eddy viscosity, which must 
be computed by a suitable turbulence model. Pr is the laminar Prandtl number, which is taken 
as 0.7 for air, Pr, is the turbulent Prandtl number, taken as 0.9, and y is the ratio of specific 
heats of the fluid. 

The high Reynolds number k-z turbulence model originally described by Launder and 
Spalding [11], can similarly be written as 
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where w , f c and g c are now given by 
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The diffusive fluxes /„ and g v are given by 
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and the source term h is given by 
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where the production term P and the term S in two dimensions are given by 
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The eddy viscosity is calculated as 
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and k also appears in the normal stresses in equation (5). The constants appearing in the 
above equations are given the standard values recommended in [11], i.e. 


C^ = 0.09 a* = 1.0 a e = 1.3 C, = 1.44 C 2 =1.92 (16) 

These equations are coupled to the governing equations for the mean flow and exhibit a similar 
structure. Therefore, a single system of equations which simultaneously governs the flow and 
turbulence quantities may be written as 
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dt dx + dy dx + dy + 
where the solution vector and the source term are now given by 
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and the flux definitions follow from equations (2),(4),(11), and (12). 
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The solution procedure consists of discretizing these equations in space on an unstruc- 
tured mesh, and then integrating the discretized equations in time until the steady-state solution 
is obtained. The basic strategy pursued in this work involves the use of a finite-element Galer- 
kin discretization technique, in conjunction with an unstructured multigrid integration technique 
to solve for the steady-state. Although all six equations of the governing system are solved 
simultaneously in the multigrid strategy, the flow equations are only loosely coupled to the tur- 
bulence equations (through the value of jx ( ), and we choose to employ somewhat different base 
grid solvers for the flow equations and the turbulence equations. 


3. SPATIAL DISCRETIZATION 

The equations governing the mean flow are discretized using a Galerkin finite-element 
approach [4], The flow variables are stored at the vertices of the triangles. The convective 
fluxes are computed at the vertices of the triangles and assumed to vary linearly over the tri- 
angular elements. For the viscous terms, the flow variables themselves are assumed to vary 
linearly over the triangular elements of the mesh, and the required velocity gradients in the 
expression for the viscous stresses are thus computed at the centers of the triangular elements. 
Additional artificial dissipation terms are required to ensure the stability of the convective 
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terms and these are constructed as a blend of a Laplacian and biharmonic operators in the con- 
served variables, designed to ensure second-order accuracy throughout the flow-field, except in 
the vicinity of a shock where first-order accuracy is recovered. For the turbulence equations, 
the diffusive terms are similarly discretized using a Galerkin finite-element approach, assuming 
linear variations of the conserved variables over the triangular elements. The velocity gradients 
in the source terms are also constructed assuming linear elements. The convective terms, how- 
ever, are constructed using first-order upwinding. Although only first-order accurate, this 
approach is employed since it helps ensure stability and positivity of the conserved variables 
throughout the integration procedure, as will be shown. Furthermore, in regions where convec- 
tion is small compared to the diffusion terms or the source terms, such as in the logarithmic 
law of the wall region, the scheme reverts to second order accuracy. In future work however, a 
second-order accurate implementation of the convective terms may be pursued. 

4. INTEGRATION SCHEME 


The discretized mean flow equations are integrated in time using an explicit five-stage 
Runge-Kutta time-stepping scheme, where the convective terms are evaluated at every stage, 
and the dissipative terms are only evaluated at the first, third and fifth stages. This scheme, 
which has previously been described [4,12], has been particularly devised to ensure rapid 
damping of high-frequency errors, and is thus well suited to drive the multigrid algorithm. 
Convergence is accelerated by the use of local time-stepping, and implicit residual averaging. 
In principle, the turbulence equations may be integrated in time using the same explicit 
scheme. However, the presence of source teims imposes a further time-step restriction. If the 
flow equations and turbulence equations are integrated in a fully coupled manner, the minimum 
local time-step from the flow and turbulence equations must be employed. In regions where the 
source terms dominate, this may lead to slow convergence. If, on the other hand, the flow 
equations and turbulence equations are integrated in an uncoupled explicit manner, the tur- 
bulence equations may significantly lag the flow equations and thus inhibit convergence to the 
steady-state solution. In order to advance the turbulence quantities at the same rate as the flow 
equations, the source terms must be treated implicitly. However, rather than simply treat the 
source terms implicitly, the system of turbulence equations is integrated in a point-implicit 
manner. Thus we rewrite the discretized turbulence equations as 


AW; 

— - /?(*,) + //(*,) (19) 

where R (w, ) represents the discretized convective and diffusive terms, which depend on the 
values of w at i and at neighboring nodes, and //(w,) represents the discretized source terms, 
which only depend on the values of w at i. The above equation is then linearized about the 
values at i which, upon solving for Aw, yields 


Aw,- = 
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The Runge-Kutta scheme described above is now replaced by a multi-stage implicit scheme, 
where the qth stage is given by 


w<«> = w<°> 


where the a q denote the Runge- 
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-Kutta coefficients for the qth stage, and At is the local time 
step. In this manner, the high-frequency damping characteristics of the original scheme are 
approximated, while the time-step restriction due to the source terms is alleviated. The precise 
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value of the local time-step At employed is one which guarantees stability as well as positivity 
of the turbulence quantities. 

5. STABILITY AND POSITIVITY CONSIDERATIONS 

One method to guarantee stability of the system is to ensure that the matrix to be inverted 
is diagonally dominant. This is not a necessary condition for stability, although it is sufficient. 
This can obviously be achieved by choosing At to be sufficiently small. However, the reason 
for employing a point-implicit approach now becomes apparent. Since the two turbulence equa- 

dR 

tions are only coupled through their source terms, the matrix is diagonal. The contribution 

from the diffusive terms is strictly negative, as well as that from the first-order upwinded con- 
vective terms. Hence, these terms, when subtracted from the diagonal of the matrix to be 
inverted, increase the diagonal dominance, and hence permit the use of a larger time-step. The 
maximum value of At is found by equating each diagonal element to its corresponding off- 
diagonal element in the coefficient matrix. The actual value employed for the time-step is 
taken as the minimum between the two values obtained by the diagonal dominance test, and 
the value determined by local stability analysis for an explicit scheme in the absence of source 
terms. 

Physically, k and e represent quantities which must remain non-negative. Thus a further 
time-step restriction is required to ensure positivity. For a simple 2x2 system, this can easily be 
derived analytically. Thus, we require that the new update to the turbulence variables be such 
that 


or, when Aw < 0 


w + Aw > aw 


0 < a < 1 


I Aw I < (1-a) w 0 < a < 1 (22) 

Substituting into equation (20), and using Cramer’s rule to evaluate the inverse of the 2x2 
matrix, we obtain two quadratic inequalities for At , i.e. one for positivity of k , and one for e. 
The time step is then limited by the smallest positive root of the two quadratic equations. 


6. MULTIGRID STRATEGY AND STEADY-STATE CONSIDERATIONS FOR THE 
k - e EQUATIONS 

A multigrid strategy is employed to accelerate the solution of the system of mean flow 
and turbulence equations to steady-state. In the context of unstructured meshes, multigrid may 
be applied by generating a sequence of non-nested coarse and fine meshes, and transferring the 
variables, residuals and corrections back and forth between the various meshes using linear 
interpolation. The patterns for interpolating between non-nested unstructured meshes are deter- 
mined in a preprocessing stage, using an efficient search algorithm. The present multigrid stra- 
tegy has previously been described in detail for the Euler and Navier-Stokes equations [4,10], 
and thus will not be repeated here. In previous multigrid applications for turbulent flows using 
algebraic models on structured and unstructured meshes [7,13], the eddy viscosities are only 
computed on the finest grid, and interpolated to the coarser meshes. Since the eddy viscosity 
represents the main coupling between the flow equations and the turbulence equations, a simi- 
lar approach has been adopted in the present context, thus ensuring a more accurate representa- 
tion of this quantity on all grid levels. However, since the eddy viscosity is only computed on 
the finest grid, it is effectively held constant throughout an entire multigrid cycle, and the 
source terms must be linearized accordingly. Making use of equations (11) and (13), the linear- 
ization of the source terms on all grids is therefore taken as 
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At this point it is worth commenting on what types of errors may be expected to be handled 
efficiently by a multigrid strategy. Multigrid is an effective device for relieving the spatial 
stiffness associated with a set of discretized equations, which is achieved by time-stepping on 
coarser grids. The turbulence equations contain spatial terms such as convection and diffusion, 
but the source terms are purely local terms. In fact, in the absence of convection and diffusion, 
the equations become completely uncoupled in space, and a properly formulated multigrid 
algorithm should yield vanishingly small corrections in such a case. Thus, it is important for 
the base grid solver to efficiently eliminate errors associated with these terms. From another 
point of view, if a purely explicit scheme were employed, a time-step restriction would arise 
from the convection, diffusion and source terms. While the first two restrictions are relaxed 
when going to coarser grids, the latter remains the same on all grid levels, effectively prevent- 
ing the use of large time-steps on coarse grids and severly limiting the overall rate of conver- 
gence. The use of a point-implicit scheme, therefore, relieves any such restrictions, and results 
in overall convergence rates similar to that achieved with the mean flow equations. 


At steady-state, the turbulence equations do not necessarily exhibit a unique solution. In 
regions where the production term \i,P vanishes, k = 0, e = 0 is an obvious solution which can 
be found by inspection of equations (10) and (13). However, the eddy viscosity, which is given 
by equation (15), becomes a ratio of two vanishing quantities, and is thus undefined. The time 
dependent turbulence equations however, are not ill-posed. On the contrary, the value of the 
constant C 2 has been carefully chosen to ensure that k, e and p, all vanish asymptotically for 
an isotropic decaying turbulence. For an isotropic turbulence, all spatial terms as well as the 
production term vanish, and equations (10) and (13) reduce to 
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Solution of this system yields the following asymptotic behavior 


k = t 
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which, for the current value of 1 < C 2 < 2 indicates that all quantities vanish for large t. 
Hence, in order to converge to the appropriate steady-state solution, it is important for any 
numerical scheme to respect the relative asymptotic time behavior of k and e throughout the 
convergence process. For the base grid solver, this is achieved by employing the maximum 
time-step for the system of two turbulence equations which ensures stability and positivity; let- 
ting i or e become negative, or taking too large time steps and subsequently limiting the 
updated values of k or e invariably leads to unrealistic values of p, in the far- field. Within the 
multigrid strategy, corrections interpolated back from the coarser grids may cause £ or e to 
become negative. Rather than limit these corrections, they are simply omitted at any point 
where positivity cannot otherwise be guaranteed. In this manner, the (point-wise) time con- 
sistency is not violated, and the overall effect is simply to lag such points by the effective 
coarse grid time step. An alternate approach would be to recompute the coarse grid corrections 
employing a smaller time step which guarantees positivity. However, due to the recursive 
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nature of multigrid, this represents a non-trivial task. 


7. BOUNDARY AND INITIAL CONDITIONS 

The derivation of the above turbulence transport equations is made under the hypothesis 
of a large Reynolds number flow. Thus, in regions close to the wall, such as in the viscous 
sublayer where molecular effects become important, these equations are not valid. In order to 
avoid integrating the turbulence equations in these region we make use of wall functions. In 
this approach, the governing equations for the flow and the turbulence are integrated up to a 
distance y =y n away from the wall. The flow in the remaining region 0 < y < y n is assumed to 
obey the law of the wall i.e. 
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At each time-step in the solution procedure of the governing equations, an estimate of the velo- 
city U at y = y n is obtained. From this, the value of the wall shear stress can be obtained by 
solving equation (26) implicitly for t/ x ( using a Newton-Raphson method). This estimate of 
the wall shear stress is then employed as a boundary condition on the momentum equation for 
the mean flow, and results in Dirichlet wall boundary conditions for k and e. In practice the 
point y = y* is very close to the wall so that it may be approximately placed on the wall, and 
the boundary conditions at y = y n may be imposed at the wall surface. For the momentum 
equation, this results in a wall slip velocity U = U(y n ). 


In the far-field, k and e are assigned freestream values at inflow boundaries, and simple 
extrapolation is employed at outflow boundaries. Initial conditions on k and e are obtained by 
imposing a level of freestream turbulence from which k is determined, and e is evaluated from 
equation (15) in order to produce a low value of freestream eddy viscosity (p, < 1). However, 
since the present formulation results in a small value of p, in all regions of the flow field 
where production is negligible, the converged solution is relatively insensitive to the initial 
values of k and e. The mean flow equations are initialized using uniform freestream flow con- 
ditions, and applying the tangential slip velocity boundary condition (as for an inviscid flow). 
Throughout the integration process, the wall shear stress obtained from equation (26), which is 
fed back into the momentum equation, retards the flow near the wall, thus creating a boundary 
layer profile. 


8. RESULTS USING WALL FUNCTIONS 

Two attached flow cases have been computed using the multigrid implementation of the 
high-Reynolds-numbcr turbulence model described above. The first case consists of transonic 
flow past an RAE 2822 airfoil. The freestream Mach number is 0.729, the incidence is 2.31 
degrees, and the Reynolds number is 6.5 million. This case (case 6) has been well documented 
both experimentally [14] and computationally [7,13,15] on structured and unstructured meshes. 
The mesh employed for this case is depicted in Figure 1. It contains 12,823 vertices and exhi- 
bits a normal spacing of 10^* chords at the airfoil surface, which positions the first point off the 
wall in the logarithmic law of the wall region. The computed Mach contours for this case are 
shown in Figure 2, while the resulting eddy viscosity distribution is given in Figure 3. A 
smooth distribution of eddy viscosity throughout the boundary-layer and wake regions, and 
vanishingly small values in the inviscid regions of flow are observed. The computed surface 
pressure distribution is compared with experimental data [14] in Figure 4, showing good 
overall agreement. The convergence rate of the system of equations is depicted in Figure 5, by 
plotting the RMS average of the density residual, and the residual of p£ throughout the flow 
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field, versus the number of multigrid cycles. As can be seen, the flow equations and turbulence 
equations converge with the same asymptotic rates. The residuals are reduced by roughly 6 
orders of magnitude over 200 cycles, yielding an overall convergence rate of 0.93 

The second case involves flow over a high-lifting four-element airfoil. This case is useful 
in demonstrating the complex geometries and resulting flow-fields which can be handled by the 
present methodology. The mesh employed is depicted in Figure 6 . It contains a total of 51,100 
vertices and a normal spacing of lx 10" 4 chords off the wall for each airfoil element. The com- 
puted Mach contours are shown in Figure 7, while the resulting eddy viscosity distribution is 
given in Figure 8 . The ease with which multiple wakes and confluent boundary layers may be 
handled by the present approach is evident from the figures. The computed surface pressure 
distribution is seen to compare favorably with experimental wind-tunnel data from [16] in Fig- 
ure 9. It should however be pointed out that such favorable agreement is in large part due to 
the attached nature of the flow. The multigrid convergence rates of the density equation and the 
k equation are depicted in Figure 10, where both equations are seen to achieve approximately 
the same asymptotic rates, decreasing by 4 orders of magnitude over 300 cycles. 


9. LOW REYNOLDS NUMBER TURBULENCE MODEL MODIFICATIONS 


While the use of the high-Reynolds-number turbulence equations in conjunction with wall 
functions is useful for a large class of wall bounded flows, it is nevertheless limited to flows 
where a logarithmic law of the wall region exists, and is thus strictly not valid for separated 
flows. An alternative approach is to modify the turbulence equations in order to account for 
low-Reynolds-number effects. Many such modifications have been proposed over the years 
with varying degrees of success [17]. One common feature of all such modifications is that 
they have proved exceedingly difficult to integrate numerically very close to the wall. The aim 
of the present work is to develop an efficient and robust technique for integrating such models, 
rather than reformulating or advocating any one model in particular. With this in mind, we 
chose to implement the simplest possible low-Reynolds-number model that has been demon- 
strated to produce good results for simple problems, with possible extensions to more complex 
models should the original version prove inadequate for more complicated flows. To this end, 
the modifications proposed by Speziale, Abid and Anderson [18] have been implemented. The 
modified turbulence equations, now given in vector form, can be written as 
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with boundary conditions at the wall given by 

k = 0 
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and employing the following values for the constants 


(27) 


( 28 ) 



- 9 - 


2 —Ke, 

Cy, = 0.09 O k = 1.36 a e = 1.36 C, = 1.44 C 2 = 1.83 [1 - ^ exp(— ) ] 
where Re, = p— — is the turbulence Reynolds number. As cun be seen, no extra source terms 

\IZ 

are introduced, and only two damping functions are required. This model is similar in form to 
the Lam-Bremhorst model [19], with the notable difference that all damping functions depend 
solely on y + . The evaluation of such functions requires the knowledge of the distance of each 
point from the closest wall. In the context of unstructured meshes, this information can be con- 
structed through the use of a generalized distance function, as outlined by Barth [20]. As with 
most low-Reynolds-number turbulence models, the current form of the model has been 
reported to be extremely stiff in near-wall regions, generally requiring the prescription of initial 
profiles in k and e in order to guarantee convergence to steady-state. Such techniques are con- 
sidered impractical for complex aerodynamic flows, and thus a more robust solution strategy 
has been pursued. The difficulties associated with the near-wall regions can be assessed by 
inspection of equations (27). When the wall boundary condition k = 0 is substituted into the e 
equation, it is seen to result in a singularity, since k appears in the denominator of this equa- 
tion. Since / 2 also vanishes at the wall, this singularity is in principle removable. However, 
the numerical integration of the e equation in its present form will only be well behaved if / 2 
and k have the same asymptotic behavior near the wall throughout the integration procedure, 
thus the need for startup profiles. The approach taken in this work is to remove the singularity 
by solving for a new variable defined as 


k — k f 2 or k = (29) 

/ 2 

Upon substituting this expression into equation (27), and using the chain rule to evaluate the 
gradient operators, one obtains the new set of equations 
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While the k equation now looks rather complicated, the new terms are only significant in the 
region / 2 < 1, and in fact, although all terms have been included, only the £V 2 / 2 term has 
been found to have a significant effect on the overall solution. At the wall, we have 

/ 2 = 0, V/ 2 = 0, V 2 / 2 > 0 

as well as 


p,=0 P = 0 S=0 «=0 

The boundary condition k= 0 implies k bounded at the wall. Since f 2 which appears as a 
denominator in the right-hand side of the k equation vanishes, we require the non vanishing 
terms in the numerator to sum to zero, thus yielding the condition 


k = 


P £ 

pV 2 / 2 


(31) 
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Upon substituting this expression into the e equation, with u = 0, p, = 0, and S = 0 one obtains 
a modified Helmoltz equation for e at the wall in the steady-state 


V p Ve 

j'* _ r 


- C 2 pV 2 / 2 e = 0 


This equation is well behaved and simple to integrate numerically, 
employed for e is taken to be 


(32) 

The boundary condition 


de 

to m 0 < 33 > 

While it is realized that this condition may not be entirely accurate at the wall [21], it is used 
at this initial stage for simplicity and may be modified in further work. 

In regions removed from the wall, the e equation remains well behaved. The k equation 

on the other hand contains the source term k V 2 / 2 . V 2 / 2 , which can be approximated as — — 

dy 2 

has the following properties 


V 2 / 2 > o y + < 3.4 (34) 

V 2 / 2 < 0 y* > 3.4 

where y + = 3.4 represents the point of inflection in the / 2 function. In regions where V 2 / 2 is 
negative or zero, the k equation is well behaved. However, V 2 / 2 large and positive represents a 
growing source teim, which can be numerically unstable. However, since the point y + = 3.4 is 
very close to the wall, and within the laminar sublayer, k can be approximated by the relation 


or 


k + = constant . y +2 


(35) 


* = £«.« ; — (36) 

n - ex P (-^-)] 2 

4.9 

which from direct simulations [22], is generally known to be valid up to y + = 10. Finally, in 
regions far away from the wall, the damping functions become unity, their derivatives all van- 
ish, and the original high-Reynolds-number equations are recovered, albeit with the new values 
of the constants advocated in [18]. Thus, in summary, the e equation given in the form (30) is 
employed throughout the entire flow-field, except at the wall, where the form (32) is used. For 
the k equation, the form given by equation (30) is employed from the far-field up to y + = 3.4, 
which is within the laminar sublayer. Below this value of y + , equation (36) is employed with 
the boundary condition for k. given by equation (31). 

The multigrid strategy previously described for the high-Reynolds number turbulence 
equations carries over in a straight-forward manner. The linearization of the source terms is 
now taken as 


dH 
d w 


2 c , Ifi V /2 
- -S + (p + — )— — 

3 Sk Pf 2 
CiC i J i J 2 P +C 2 ^ 


fz 

-}c,S-2C,f 


(37) 


where the production term in the e equation has been simplified by the definition of p, in equa- 
tions (27), in order to remove k from the denominator. The damping functions are evaluated 
only on the finest grid, and interpolated up to the coarser grids, thus affording a more con- 
sistent representation of the equations on all grid levels. A full multigrid strategy is employed, 


where grid sequencing is used to provide an initial solution for the fine grid. In general, it has 
been found advantageous to use the high-Reynolds-number model with wall functions on 



- 11 - 


coarse grids, and the low-Reynolds-number model on the finest grid when grid sequencing, 
thus rapidly setting up appropriate levels of eddy viscosity on the finest grid. 

10. RESULTS 

The present implementation of the low-Reynolds-number turbulence model has been 
employed to compute the turbulent boundary layer over a flat plate, the transonic flow over an 
RAE 2822 airfoil, and the transonic flow over a two-element airfoil. 

The mesh employed to compute the flat plate boundary layer case is depicted in Figure 

11. It contains 24 points ahead of the plate, 48 points along the plate in the streamwise direc- 
tion, and 80 points in the direction normal to the plate. The freestream Mach number is 0.3, 
and the Reynolds number of the flow, based on the length of the plate is 10 million. The first 
point normal to the plate is located at a distance of 2x10 plate lengths, which lies in the 
region y + <l. The resulting velocity profiles are plotted in Figures 12 and 13, both in physical 
coordinates, and logarithmic wall coordinates, and compared with the well known l/7th power 
law distribution, and logarithmic law of the wall profile. The computed skin friction is plotted 
in Figure 14, versus the experimental data taken from [23], The resulting distributions of k 
and e are shown in Figures 15 and 16. The well known peaks of k and e are observed, and a 
non-zero value of e at the wall is obtained. These distributions are however slightly different 
from those obtained previously with the same model [18], and may be attributed either to the 
different boundary condition, or to the near-wall grid resolution. The overall flow quantities are 
nevertheless well predicted, as shown in Figures 12 through 14. 

The transonic flow case over the RAE 2822 airfoil presented in the previous section has 
been recomputed with the low-Reynolds-number turbulence model (Mach = 0.729, Incidence = 
2.31 degrees, Re = 6.5 million). The mesh employed is similar to that shown in Figure 1, 
except that the normal spacing at the wall is now reduced to lxlO -5 chords, which results in 
the first mesh point off the wall in the region 1 < y + < 3 over the entire surface of the airfoil. 
The computed Mach contours and eddy viscosity contours are similar to those depicted in Fig- 
ures 2 and 3, except in the near-wall regions, where both quantities vanish rapidly. The com- 
puted surface pressure and skin-friction distributions are compared with experimental data in 
Figures 17 and 18. The computed lift is slightly lower than that predicted with the wall func- 
tions and that previously obtained using an algebraic model [7]. At present, it is not clear 
whether this is due to the actual model formulation, or is associated with the present imple- 
mentation (artificial dissipation, grid resolution). However, the differences are rather small and 
the skin friction appears to be well predicted. The convergence of the density equation and the 
two turbulence equations is depicted in Figure 19, where the residuals are plotted versus the 
number of multigrid cycles on the finest grid. The flowfield and turbulence equations are all 
initialized with uniform freestream values, and 25 cycles were performed on the previous 
coarser grid using wall functions, prior to initializing the solution procedure on the finest grid. 
Initializing the calculation with freestream values for all equations on the finest grid has also 
been employed with little degradation in convergence. From Figure 19, all equations are seen 
to converge at approximately the same rate, resulting in a residual reduction of 4 to 5 orders of 
magnitude over 300 multigrid cycles. 

The final case involves the transonic flow over a two-element airfoil. This case illustrates 
the ease with which complex geometries and flows with multiple viscous layers may be han- 
dled by the present methodology. The mesh employed is depicted in Figure 20. It contains a 
total of 28,871 vertices, with a normal spacing of 2xl0 -5 chords off the wall for each airfoil 
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element. The freestream Mach number is 0.5, the incidence is 7.5 degrees, and the Reynolds 
number is 4.5 million. The computed Mach contours and eddy viscosity contours are depicted 
in Figures 21 and 22. At these conditions, the flow is supercritical and a shock is formed on 
the upper surface of the slat. A small region of separated flow occurs behind the shock, as can 
be seen from the skin friction plot of Figure 23. This region of separation has previously been 
reported in prior calculations using an algebraic turbulence model [7], The computed surface 
pressure distribution is seen to compare favorably with experimental wind-tunnel data [24], in 
Figure 24. The convergence rate for this case is depicted in Figure 25, where the residuals of 
the density equation and the two turbulence equations are reduced by approximately 3 to 4 ord- 
ers of magnitude over 300 cycles on the finest grid. 

11. CONCLUSION 

A multigrid strategy for solving the steady-state high and low-Reynolds number k - e tur- 
bulence equations has been formulated and implemented on unstructured meshes. A variety of 
aerodynamic flows have been computed, consistently demonstrating similar convergence rates 
for the turbulence and flow equations. Initialization of all flow and turbulence quantities may 
be performed using uniform freestream values. At present, the evaluation of the turbulence 
terms requires a significant fraction of the overall time within a single time-step. For example, 
the RAE 2822 supercritical airfoil flow case with the low-Reynolds-number turbulence model 
requires roughly 2.5 seconds per multigrid cycle on a single processor of the CRAY-YMP 
supercomputer, which is almost 75% higher than that required by the algebraic model reported 
previously [7], However, it is estimated that this can be substantially reduced by assembling 
the turbulence and flow residuals simultaneously within a single loop. Given the demonstrated 
convergence rates, the two-equation turbulence model should be competitive in terms of com- 
puter resources with algebraic models, while providing much greater flexibility in dealing with 
complex geometries and flow-fields. 

Future work should involve a more thorough investigation of the various two-equation 
turbulence models and their ability in predicting complex aerodynamic flows, including flows 
with massive separation. 
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Figure 1 

Unstructured Mesh Employed for Computing Flow Over an RAE 2822 Airfoil 
(Number of Vertices = 12,823) 



figure 2 

Computed Mach Contours for Turbulent Row over RAE 2822 Airfoil 
Using Wall Functions (Mach = 0.729, Re = 6.5 million. Incidence = 2.31 degrees) 
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Figure 3 

Computed Eddy Viscosity Contours for Turbulent Row over RAE 2822 Airfoil 
Using Wall Functions ( Mach = 0.729, Re = 6.5 million, Incidence = 2.31 degrees ) 
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Figure 4 

Comparison of Computed Surface Pressure using Wall Functions with Experimental 
Measurements for Flow over an RAE 2822 Airfoil 
(Mach = 0.729, Re = 6.5 million, Incidence = 2.31 degrees) 
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Number of Cycles 


Figure 5 

Convergence Rate of Density Equation and K equation Using Wall Functions 
versus the Number of Multigrid Cycles for Flow over an RAE 2822 Airfoil 




Figure 6 

Unstructured Mesh Employed for Computing Flow Over a Four-Element Airfoil 
(Number of Vertices = 51,100) 


Figure 7 

Computed Mach Contours Using Wall Functions for Flow over a Four-Element Airfoil 
(Mach = 0.2, Re = 2.83 million. Incidence = 8.18 degrees) 





Figure 8 

Computed Eddy Viscosity Contours Using Wall Functions for How over a Four-Element Airfoil 
(Mach = 0.2, Re = 2.83 million, Incidence = 8.18 degrees) 
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Figure 9 

Comparison of Computed Surface Pressure using Wall Functions with Experimental 
Wind-Tunnel Data for Flow over a Four-Element Airfoil 
(Mach = 0.2, Re = 2.83 million, Incidence = 8.18 degrees) 
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Figure 10 

gence Rate of Density Equation and K Equation Using 
inctions for Flow over a Four-Element Airfoil 
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Figure 11 

Triangular Mesh Employed for Flat Plate Boundary Layer Calculation 
(Number of Vertices = 5913, 10:1 Magnification in Y-direction) 
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Figure 14 — ' 

Computed Skin Friction Distribution for Hat Plate Boundary Layer 
Versus Experimental Data from [23] 
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Figure 16 

Computed Near Wall Distribution of e + for Flat Plate Boundary Layer 
(Mach = 0.3, Re, = 5.3 million) 
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Figure 17 

Computed Surface Pressure Distribution Using Low-Reynolds Number Modification 
for Turbulence Equations Versus Experimental Data for Flow past RAE 2822 Airfoil 
(Mach = 0.729, Re = 6.5 million, Incidence = 2.31 degrees) 



Figure 18 

Computed Skin-Friction Distribution Using Low-Reynolds Number Modification 
for Turbulence Equations Versus Experimental Data for Flow past RAE 2822 Airfoil 
(Mach = 0.729, Re = 6.5 million, Incidence = 2.31 degrees) 
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Figure 20 

Global View of Coarse Unstructured Mesh and Gose-Up View of Fine 
Unstructured Mesh Employed for Computing Flow Past a Two-EIemcnt Airfoil 
(Coarse Mesh Points = 7272, Fine Mesh Points = 28871) 






Figure 21 

Computed Mach Contours Using Low-Reynolds Number Modification for Turbulence 
Equations for Supercritical Flow over a Two-Element Airfoil 
(Mach = 0.5, Re = 4.5 million, Incidence = 7 .5 degrees) 
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Figure 22 

Computed Eddy Viscosity Contours Using Low-Reynolds Number Modification 
for Turbulence Equations for Supercritical Flow over a Two-Element Airfoil 
(Mach = 0.5, Re = 4.5 million. Incidence = 7.5 degrees) 
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Figure 24 

Computed Surface Pressure Distribution Using Low-Reynolds Number Modification 
for Turbulence Equations for Supercritical Flow over a Two-Element Airfoil 
Versus Experimental Wind Tunnel Data 
(Mach = 0.5, Re = 4.5 million, Incidence = 7.5 degrees) 
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Figure 25 

Multigrid Convergence Rate of the Density Equation and the Two Turbulence 
Equations Using Low-Reynolds Number Modifications for Flow Over 
Two-Element Airfoil (Mach = 0.5, Re = 4.5 million, Incidence = 7.5 degrees) 
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